Cite this Article
Liang Tu-Lu, Shao Wei, Wei Xiao-Kun, Liang Mu-Sheng. Hybrid sub-gridding ADE–FDTD method of modeling periodic metallic nanoparticle arrays. Chinese Physics B, 2018, 27(10): 100204  
Permissions
Hybrid sub-gridding ADE–FDTD method of modeling periodic metallic nanoparticle arrays
Liang Tu-Lu, Shao Wei †, Wei Xiao-Kun, Liang Mu-Sheng
School of Physics, University of Electronic Science and Technology of China, Chengdu 610054, China

 

† Corresponding author. E-mail: weishao@uestc.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 61471105 and 61331007).

Abstract

In this paper, a modified sub-gridding scheme that hybridizes the conventional finite-difference time-domain (FDTD) method and the unconditionally stable locally one-dimensional (LOD) FDTD is developed for analyzing the periodic metallic nanoparticle arrays. The dispersion of the metal, caused by the evanescent wave propagating along the metal-dielectric interface, is expressed by the Drude model and solved with a generalized auxiliary differential equation (ADE) technique. In the sub-gridding scheme, the ADE–FDTD is applied to the global coarse grids while the ADE–LOD–FDTD is applied to the local fine grids. The time step sizes in the fine-grid region and coarse-grid region can be synchronized, and thus obviating the temporal interpolation of the fields in the time-marching process. Numerical examples about extraordinary optical transmission through the periodic metallic nanoparticle array are provided to show the accuracy and efficiency of the proposed method.

PACS: 02.70.Bf;02.60.Cb;92.60.Ta
Keyword:locally one-dimensional finite-difference time-domain;metallic nanoparticle;sub-gridding;surface plasmon polaritons
1. Introduction

Since the extraordinary optical transmission (EOT) phenomenon through a two-dimensional (2D) metallic hole was first reported by Ebbesen in 1998,[1] a great deal of attention has been paid to the investigation of surface plasmon polaritons (SPPs) excited in sub-wavelength metallic structures.[2–9] Generally speaking, the transmission resonances in two cases, namely the coupled SPP resonant mode and the Fabry–Pérot (FP)-like resonant mode, have been involved in the interpretation of the EOT phenomenon. For the analysis of periodic metallic nanoparticle arrays, the frequency-dependent finite-difference time-domain (FDTD) method has been widely used.[10–13] Because the SPPs are highly localized along the metal-dielectric interface, fine spatial grids are required to obtain sufficient accuracy. Thus, an extremely small time step constrained by the Courant–Friedrich–Levy (CFL) stability condition results in a long CPU time.[14] Several unconditionally stable FDTD methods,such as the alternating-direction implicit (ADI) FDTD,[15,16]the split-step (SS) FDTD,[17–19] the locally one-dimensional (LOD) FDTD,[20–22] the Crank–Nicolson (CN) FDTD method,[23,24] and the Newmark-Beta FDTD method,[25] have been developed to eliminate the CFL condition in the past two decades. For the ADIFDTD and LODFDTD methods, which only need to solve a tri-diagonal matrix equation with low computational complexity, the latter requires fewer arithmetic operations and less execution time than the former.[26]

An unconditionally stable FDTD method can be utilized to simulate periodic metallic nanoparticle arrays to improve the efficiency. However, its time step cannot be chosen to be too large due to the numerical dispersion error.[27,28] More importantly, the fine grid division along the metal-dielectric interface leads to a large memory requirement. A sub-gridding scheme (sub-gridding means the fine grids inside the coarse ones) for FDTD only introduces fine grids into the regions that need high resolution.[29] A more efficient scheme, based on the use of ADI-FDTD in the local fine-grid region,[30] synchronizes the time marching of the fine-grid region with that of the coarse-grid region. Thus, only spatial interpolation of the fields on the spatial sub-gridding interface is necessary.

In this work, we propose a modified sub-gridding scheme that hybridizes FDTD and LOD-FDTD for the EOT analysis of periodic metallic nanoparticle arrays. After the spatial domain is divided into the local fine-grid region and global coarse-grid region, the LOD-FDTD is applied to the fine-grid region while FDTD is used for the rest of the coarse-grid region. An easy-to-accomplish scheme of spatial interpolation is proposed at the interface between coarse-grid and fine-grid regions, and its numerical stability is tested. Meanwhile, the time step for fine-grid regions can be taken to be the same as that for coarse-grid regions due to the unconditional stability of LOD-FDTD. The dispersion of the metal is expressed by the Drude model, which is solved with the auxiliary differential equation (ADE) technique[31] to establish the relationship between the electric field intensity and conductive electric current in metallic nanoparticles. With the hybrid sub-gridding ADE–FDTD method, the transmission resonances of periodic metallic nanoparticle arrays with and without the defect structure are investigated from the nearultraviolet region to the nearinfrared region. The numerical examples verify the accuracy and the effectiveness of the proposed method.

2. Hybrid sub-gridding ADE–FDTD method

The time dependence of e−iωt is assumed. According to the Drude model, the relative permittivity (in frequency domain) of the dispersive metal is given by

where ω is the angular frequency of the impinging light, ωp is the plasma frequency in the metal, ε∞ is the relative permittivity at infinite frequency in the metal, and γ is the damping constant in the metal.

For simplicity, a 2D TM wave including the Ex, Ez, and Hy components is considered. With reference to Ref. [31], the Drude model is handled by using the ADE technique, and the following time-domain governing equations are obtained, including the Maxwell equations and auxiliary differential equations:

where εeff = ε for free space and εeff = ε ε∞ for the metal region.

The 2D solution domain is divided into the local fine-grid region and the global coarse-grid region, as shown in Fig. 1. The ADE–LOD–FDTD is then applied to the fine-grid region while ADE–FDTD is for the remaining coarse-grid region.

Fig. 1. Diagram of a sub-gridding scheme for a grid refinement ratio of three.

Because of unconditional stability of LOD-FDTD, a large time step can be taken for the fine grids. In other words, the same time step can now be used in both the coarse-grid region and the fine-grid region. One of the immediate benefits is that an interpolation scheme needs to be developed in the spatial domain instead of the time domain; the other benefit is that such a scheme reduces the memory requirement since the relatively memory-expensive LOD-FDTD is only applied to the fine-grid region.

2.1. Coarse-grid region

The ADE–FDTD method is applied to the coarse-grid region, and we have (taking Ez and Jx for example, where the fields and currents in coarse grids are shown in capital letters)

where ΔX is the mesh size in the coarse-grid region.

Since the infinite periodicity assumed in this paper lies in the x direction, the periodic boundary condition (PBC) is imposed on the top and bottom boundaries, as shown in Fig. 1.

2.2. Fine-grid region

In the ADE–LOD–FDTD fine-grid region, the well-known Yee’s finite difference grid is used. The electromagnetic-field components are arranged on the cells in the same way as that using the ADE–FDTD method.

With the LOD scheme,[20] we obtain the ADE–LOD–FDTD formalism. In the first time step (n + 1/2), we have (taking ez and jx for example, where the fields and currents in fine grids are shown in lower case letters)

where Δx is the mesh size of the fine-grid region. The time step is taken uniformly to be the same in both coarse-grid and fine-grid regions.

The equations of the second time step (n + 1) can be obtained in a similar way. The above two steps are performed alternatively in the unconditionally stable time-marching procedure of 2D ADE–LOD–FDTD. A more detailed process has been introduced in Ref. [32].

2.3. Implementation of sub-gridding

The interpolation between field samples of coarse and fine grids affects the stability and accuracy of the hybrid sub-gridding ADE–FDTD. In Ref. [33] a sub-gridding scheme for the 2D case was proposed, which takes on low reflections from the spatial sub-gridding interface.

Figure 2 shows the spatial sub-gridding interface between the coarse-grid region and the fine-grid region. With field samples from Fig. 2, the electric and magnetic field components by using linear interpolation and the collocated field samples are given as follows:[33]

Fig. 2. Arrangement of electromagnetic field components at the spatial sub-gridding interface.

In order to further improve the accuracy and stability of the sub-gridding ADE–FDTD method, a modified formula for evaluating the magnetic field of coarse grids in the sub-gridding region is proposed below.

Thus, the electric fields of fine grids in the spatial sub-gridding interface and its adjacent co-directional electric fields are calculated by using Eq. (11), and magnetic fields of coarse grids in the sub-gridding region are calculated with Eq. (13).

To test the stability of the proposed sub-gridding scheme, we calculate a 2D square resonator of 255 nm × 255 nm with Δfine = 5 nm and Δcoarse = 3Δfine nm. The results shown in Fig. 3 indicate that our sub-gridding scheme is still stable after 100 thousand time-marching steps. The scheme of spatial interpolation introduced in Ref. [33] shows unstable behavior after 5000 time-marching steps.

Fig. 3. (color online) Numerical stability verification of hybrid sub-gridding ADE–FDTD: (a) scheme of spatial interpolation in this paper, (b) scheme of spatial interpolation in Ref. [28].

The basic flowchart of the hybrid sub-gridding ADE–FDTD is shown in Fig. 4.

Fig. 4. Flowchart of the hybrid sub-gridding ADE–FDTD method.
3. Numerical results

With the proposed hybrid sub-gridding ADE–FDTD method, the transmission resonances of periodic metallic nanoparticle arrays with and without defect structures are investigated from the nearultraviolet region to the near infrared region in this section.

A normally incident TM-polarized plane wave illuminates a periodic metallic nanoparticle array from the left side, as shown in Fig. 5. A Gaussian pulse is used as the source excitation, which can be written as

where the maximum frequency fmax = 1000 THz, τ = 1/(2fmax), and t0 = 3τ. The frequency characteristics of the transmission are calculated with the discrete Fourier transform (DET) of time-domain response in the observation plane. The metal investigated in this paper is gold and the considered wavelength range is from 300 nm to 1500 nm. Here we choose ε∞ = 1, ωp = 1.37 × 1016 rad/s, and γ = 4.08 × 1013 rad/s as a set of parameters to fit the experimental relative permittivity data from the nearultraviolet region to the nearinfrared region.[34]

Fig. 5. Schematic diagram of the periodic metallic nanoparticle array without defect structure.

The computational region is truncated by Berenger’s perfectly matched layer (PML) on the left and right boundaries. The top and bottom boundaries are treated by the PBCs due to the periodicity of the structure.

The computational domain is discretized with square Yee grids. For the fine-grid region simulated with ADE–LOD–FDTD, we choose Δfine = 2.5 nm. For the coarse-grid region with ADE–FDTD, we choose Δcoarse = 3Δfine. The time step size is set to be ΔtFDTD = ΔtLOD = Δcoarse/2c (c is the velocity of light in a vacuum) for both ADE–FDTD and ADE–LOD–FDTD. Due to the numerical dispersion error,[28] a large ratio of Δcoarse to Δfine will degenerate the accuracy of the hybrid sub-gridding ADE–FDTD method.

3.1. Accuracy verification of hybrid sub-gridding ADE–FDTD method

First, we investigate the light transmission through the periodic metallic nanoparticle array shown in Fig. 5. In this figure, r, d, and T correspond to the radius of the nanoparticle, the space between two adjacent nanoparticles, and the period length of the array, respectively.

Figure 6 shows the transmittance values calculated by the hybrid sub-gridding ADE–FDTD method and the fine-grid ADE–LOD–FDTD method, where r = 75 nm, d = 45 nm, and T = 210 nm. It can be seen from the figure that the results from the two methods are in good agreement. Table 1 shows the comparison of computational efforts between the two methods. With the sub-gridding scheme, the CPU time of the proposed method can be much less than that of the fine-grid ADE–LOD–FDTD. Since the relatively memory-expensive ADE–LOD–FDTD is only applied to the fine-grid region, the memory requirement of the hybrid sub-gridding ADE–FDTD method is about half that of the fine-grid ADE–LOD–FDTD. All calculations in this paper are performed on an Intel Core (TM) i7 3.60-GHz computer with 8-GB RAM.

Fig. 6. (color online) Plots of transmittance versus wavelength from two methods.
Table 1.

Comparison of computational efforts between two methods.

.
3.2. Effects of radius of nanoparticles

With the proposed hybrid sub-gridding ADE–FDTD method, the results of the transmission spectrum for nanoparticle radius variation at the fixed parameters of d = 45 nm and T = 210 nm are depicted in Fig. 7. The radius of the nanoparticle is changed from 60 nm to 90 nm. As the radius of the nanoparticle increases, the resonance peak has an obvious red shift in the transmission spectrum, and the photonic band gap in the periodic metallic nanoparticle array is broadened to high frequency.

Fig. 7. (color online) Plots of transmittance versus wavelength for different nanoparticle radii.

It can be seen that the resonance peak is very sensitive to the change of the radius of the nanoparticle. Therefore, we can choose the appropriate radius of the nanoparticle so that the resonance peak appears at a desired wavelength range.

3.3. Effect of spacing between nanoparticles

The effect of the nanoparticle spacing variation is illustrated in Fig. 8 when r = 75 nm and T = 210 nm. The spacing between two adjacent nanoparticles varies from 22.5 nm to 75 nm.

Fig. 8. (color online) Plots of transmittance versus wavelength for different spacing between nanoparticles.

As the spacing between two adjacent nanoparticles increases, the decrease of the energy to drive the resonance results in an obvious red shift of the resonance peak in the transmission spectrum.

3.4. Effects of period length of array

Figure 9 depicts the results of the transmission spectrum for periodic length variation at the fixed parameters of r = 75 nm and d = 45 nm. The period length of the array is changed from 180 nm to 240 nm.

Fig. 9. (color online) Plots of transmittance versus wavelength for different period lengths of array.

As the period length of the array increases, the resonance peak has an obvious blue shift in the transmission spectrum, and the photonic band gap of the periodic metallic nanoparticle array is compressed towards low frequencies.

3.5. Effect of defect structure

By adding point defects or line defects into a periodic nanoparticle array to change its optical properties, the structure can be applied to the manufacture of various optical devices such as filters, optical waveguides, light-emitting diodes, etc. In this part, we discuss the effect of defect nanoparticle size on transmission of a periodic metallic nanoparticle array.

Figure 10 shows a schematic view of one period of the 2D periodic metallic nanoparticle array with defect structures. The defect nanoparticle is placed in the center of the array with a radius of R, and the other four identical nanoparticles each have a radius of r. With the proposed hybrid sub-gridding ADE–FDTD method, the results of the transmission spectrum for defect nanoparticle radius variation are shown in Fig. 11. The radius of the defect nanoparticle changes from 60 nm to 90 nm, while the following parameters are fixed: d = 45 nm, r = 75 nm, and T = 210 nm.

Fig. 10. (color online) Schematic diagram of periodic metallic nanoparticle array with defect structure.
Fig. 11. (color online) Plots of transmittance versus wavelength of periodic metallic nanoparticle array with defect structures.

As the radius of the defect nanoparticle increases, the photonic band gap shifts toward the longwavelength region. The red shift phenomenon is obvious, and the photonic band gap is broadened. When no defect is in the array where the radius of the defect nanoparticle is 75 nm, the transmission peak is the largest. This phenomenon is caused by the interaction between the defect nanoparticle and its adjacent ones. No matter whether the radius of the defect nanoparticle increases or decreases, the amplitude of the transmission peak decreases.

The introduction of the defect nanoparticle into the periodic array broadens the bandwidth of the photonic band gap. This shows that the resonant wavelength of the transmission peak is tunable, and these human-controlled factors can greatly expand the application scope of periodic metallic nanoparticle arrays, such as wavelength-selective devices, small-size nonlinear optical sensors, and filters.

3.6. Effect of bending waveguide

With the proposed hybrid sub-gridding ADE–FDTD method, in this part, we discuss the effect of a bending waveguide composed of periodic metallic nanoparticles. The radius of the metallic nanoparticle is r = 75 nm, and the space between two adjacent nanoparticles is d = 45 nm.

Figure 12 shows the transmission curve obtained from the bending waveguide formed by two sharp corners with a 90° angle. From the figure, one can see that most of the guided modes are transmitted except that a transmission dips at a wavelength of 420 nm. Figure 12 shows a numerical illustration of the propagation of the wave except 420 nm passing through the bending waveguide in which the incident wave propagates along the first straight waveguide, couples successfully with the perpendicular one, then reaches to the last horizontal one.

Fig. 12. (color online) Plots of transmittance versus wavelength of a bending waveguide composed of periodic metallic nanoparticles.
4. Conclusions

In this paper, an efficient modified hybrid sub-gridding ADE–FDTD method based on a combination of FDTD and LOD–FDTD has been developed for the EOT analysis of periodic metallic nanoparticle arrays. The LOD-FDTD is applied to the local fine-grid region, and the FDTD is applied to the global coarse-grid region. An important advantage of this method is that the LOD–FDTD fine-grid region can be synchronized with the time-marching step of the FDTD coarse-grid region. Thus, an interpolation scheme needs to be developed in the spatial domain instead of the time domain, and such a scheme reduces the memory requirement and CPU time. A modified interpolation scheme that is easy to implement is proposed for the spatial sub-gridding interface, and its good numerical stability is proved.

We apply the hybrid sub-gridding ADE–FDTD method to the computational domain that includes general dispersive media and curved boundaries to conduct the periodic metallic nanoparticle array studies. Several periodic metallic nanoparticle arrays with and without the defect structure and a bending waveguide are also discussed from the nearultraviolet region to the nearinfrared one. The numerical results from the hybrid sub-gridding ADE–FDTD agree well with those from fine-grid ADE–LOD–FDTD, and its computing time is greatly reduced. It is found that some distinct resonance peaks appear in the photonic band gap and the resonance frequency depends strongly on the radius of the nanoparticle, the spacing between two adjacent nanoparticles, and the periodic length of the array. It is also shown that the transmittance band gap is broadened and the red shift of the transmission resonance peak occurs by introducing the defect structure.

Reference
[1] Ebbesen T W Lezec H J Ghaemi H F Thio T Wolff P A 1998 Nature 391 667
[2] Porto J A Garcia-Vidal F J Pendry J B 1999 Phys. Rev. Lett. 83 2845
[3] Min C J Wang P Chen C C Deng Y Lu Y H Ming H Ning T Y Zhou Y L Yang G Z 2008 Opt. Lett. 33 869
[4] Khanikaev A B Mousavi S H Shvets G Kivshar Y S 2010 Phys. Rev. Lett. 105 126804
[5] Moreno E Martin-Moreno L Garcia-Vidal F J 2006 J. Opt. A Pure Appl. Opt. 8 S94
[6] Martin-Moreno L Garcia-Vidal F J Lezec H J Pellerin K M Thio T Pendry J B Ebbesen T W 2001 Phys. Rev. Lett. 86 1114
[7] Liu H T Lalanne P 2008 Nature 452 728
[8] Genet C Ebbesen T W 2007 Nature 445 39
[9] Gordon R Sinton D Kavanagh K L Brolo A G 2008 Acc. Chem. Res. 41 1049
[10] Yang J L Li R P Han J H Huang M J 2016 Chin. Phys. 25 083301
[11] Li Q B Wu R X Yang Y Sun H L 2013 Chin. Phys. Lett. 30 074208
[12] Zhuan S X Ma X K 2012 Acta Phys. Sin. 61 110206 in Chinese
[13] Li L Q Shi Y X Wang F Wei B 2012 Acta Phys. Sin. 61 125201 in Chinese
[14] Taflove A Hagness S C 2000 Computational Electrodynamics: The Finite-Difference Time-Domain Method Norwood, MA Artech House
[15] Namiki T 1999 IEEE Trans. Microw. Theory Tech. 47 2003
[16] Zheng F H Chen Z Z Zhang J Z 1999 IEEE Microw. Guided Wave Lett. 9 441
[17] Lee J Fornberg B 2003 J. Comput. Appl. Math. 158 485
[18] Lee J Fornberg B 2004 J. Comput. Appl. Math. 166 497
[19] Gao L Zhang B Liang D 2007 J. Comput. Appl. Math. 205 207
[20] Shibayama J Muraki M Yamauchi J Nakano H 2005 Electron. Lett. 41 1046
[21] Nascimento V E D Borges B H V Teixeira F L 2006 IEEE Microw. Wireless Compon. Lett. 16 398
[22] Wei X K Shao W Shi S B Zhang Y Wang B Z 2015 Chin. Phys. 24 070203
[23] Sun G Trueman C W 2004 IEEE Trans. Anten. Propag. 52 2963
[24] Sun G Truneman C W 2006 IEEE Trans. Microw. Theory Tech. 54 2275
[25] Shi S B Shao W Wei X K Yang X S Wang B Z 2016 IEEE Microw. Theory Tech. 64 4082
[26] Ahmed I Chua E K Li E P Chen Z Z 2008 IEEE Trans. Anten. Propag. 56 3596
[27] Zheng F Chen Z 2001 IEEE Trans. Microw. Theory Tech. 49 1006
[28] Ahmed I Chun E K Li E P 2010 IEEE Trans. Anten. Propag. 58 3983
[29] Chevalier M W Luebbers R J Cable V P 1997 IEEE Trans. Anten. Propag. 45 411
[30] Wang B Z Wang Y Yu W Mittra R 2001 IEEE Trans. Adv. Packag. 24 528
[31] Gray S K Kupka T 2003 Phys. Rev. 68 045415
[32] Liang T L Shao W Shi S B Ou H 2016 IEEE Photon. J. 8 7804710
[33] Kulas L Mrozowski M 2008 Proc. 38th European Microwave Conf. 658 10.1109/EUMC.2008.4751538
[34] Palik E D 1982 Handbook of optical constants in solids Academic